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Abstract 

We propose a computational method to solve optimal swimming problems, 
based on the boundary integral formulation of the hydrodynamic interaction 
between swimmer and surrounding fluid and direct constrained minimization 
of the energy consumed by the swimmer. 

We apply our method to axisymmetric model examples. We consider a 
classical model swimmer (the three-sphere swimmer of Golestanian et al.) as 
well as a novel axisymmetric swimmer inspired by the observation of biolog- 
ical micro-organisms. 
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1 Swimming at low Reynolds Numbers 



Swimming is the ability to advance within a fluid by performing cyclic shape 
changes, in the absence of external propulsive forces. One of the main diffi- 
culties of swimming at small scales is given by the time-reversal property of 
Stokes flows, which describe hydrodynamics at low Reynolds numbers. 

The physical implication of this mathematical property is that simple 
swimming strategies used in nature at larger scales, for example by scal- 
lops, where the same shape change is executed forward and backward at 
different velocities to achieve propulsion in one direction, do not work for 
micro-, or nano-swimmers. At this scale, swimmers have to undergo non- 
reciprocal deformations in order to achieve propulsion. Also known as the 
Scallop-Theorem, this fact is discussed by Purcell in Ref. PQ, where a simple 
mechanical device that can indeed swim at low Reynolds numbers, the three- 
link swimmer, is proposed. A mathematical statement of the scallop theorem 
and its proof can be found in Ref. J2J. More recently several other simple 
swimmers have been presented, for example, in Golestanian and Najafi j3] 
and Avron et al. [4 J 

A mathematical approach to the problem of finding an optimal stroke 
has been proposed by Alouges et al. in Ref. [3], where it is shown how to 
formulate and solve numerically the problem of finding optimal strokes for 
low Reynolds number swimmers by focusing on the three-sphere swimmer 
of Najafi and Golestanian [3j (a simple, yet representative example). The 
analysis carried out in Ref. [5J shows how to address quantitatively swimming 
as the problem of controlling shape in order to produce a net displacement 
at the end of one stroke. By casting it in the language of control theory, 
the problem of swimming reduces to the controllability of the system, and 
the search of optimal strokes to an optimal control problem leading to the 
computation of suitable sub-Riemannian geodesies. 

The use of numerical algorithms to find optimal strokes can lead to dra- 
matic improvements. For the three-sphere swimmer, one can achieve an in- 
crease of efficiency exceeding 300% with respect to more naive proposals. j5] 
The simplicity of three-sphere swimmer, which is a system with two degrees 
of freedom, enables one to carry out the analysis in full detail. The study of 
biologically relevant swimmers, however, requires more abstract mathemat- 
ical tools and more efficient numerical algorithms. The aim of this paper is 
to introduce further numerical tools to overcome some of the computational 
limits of the numerical method used in Ref. [5], opening up the possibility 
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to treat more and more complex swimmers such as the new model swimmer 
described in section |4j which is inspired by biological micro-swimmers like 
unicellular eukaryotic algae. 

For simplicity, we focus on the special case of axisymmetric swimmers, 
which provides an interesting balance between complexity and generality of 
the attainable results. A complete theory to analyze axisymmetric swimmers 
whose shape depends on finitely many parameters and a general method to 
determine strokes of maximal efficiency is presented in Ref. [6]. The main 
feature qualifying this approach is the possibility of resolving hydrodynamic 
interactions arising from the swimmer motion in their full complexity, with- 
out being confined to asymptotic regimes (dilute limits for assemblies of 
spheres, slender body hydrodynamics for one-dimensional objects) in which 
explicit formulas for hydrodynamic forces are available. 

This work follows a similar line of thought. The main novelties intro- 
duced here are the use of a boundary integral formulation for the solution 
of the axisymmetric Stokes flow induced by the motion of the swimmer, and 
the resolution of the optimal control problem via direct constrained mini- 
mization of the functional giving the energy spent by the swimmer. This 
should be contrasted with the resolution of the hydrodynamic interactions 
via the finite element method that was used in Ref. [5] and Ref. [6] , and the 
resolution of the Euler-Lagrange equation associated with the energy min- 
imization problem via an ODE scheme coupled with a shooting method to 
enforce periodicity of the swimming stroke. 

The accuracy that can be obtained with the new method is a few or- 
ders of magnitude better than with the previous one, due to the possibility 
of resolving the swimmer details at a much smaller scale at a comparable 
computational cost. In addition to the improved accuracy granted by the 
use of boundary integrals, the direct solution of the constrained minimiza- 
tion problem enlarges the class of optimization problems we can consider. In 
particular, it enables us to consider also inequality constraints on the shape 
variables. This is often necessary in order to restrict the set of admissible 
shape changes to those that are also biologically plausible. 

The rest of this paper is organized as follows. In section [2] we derive 
the equations characterizing optimal strokes, while in section [3] we propose a 
strategy for their numerical solution which is alternative to the ones presented 
in Ref. [5] and Ref. [6] . We validate our method by comparing with previous 
results on the three-sphere swimmer, and by calculating optimal strokes for 
a new model axisymmetric swimmer in section |4} 
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2 Governing Equations 



We consider an axisymmetric micro-organism or micro-robot swimming at 
low Reynolds numbers in a three dimensional fluid at rest at infinity. No 
interaction with other obstacles is considered. 

The swimmer occupies at time t the (unknown) portion of space denoted 
by Q(t) C 9? 3 , with boundary dQ(t) = T(t), while the surrounding region 
3? 3 \f2(t) is occupied by the fluid. 

We are interested in finding the "optimal stroke" of period T, directed 
along the axis of symmetry of the swimmer, and that will take the swimmer 
from the position O to the position fi(T) = ^o + ce, where c/T is the average 
swimming velocity. 

According to whether we want to optimize performance or energy con- 
sumption, different optimality criteria can be adopted. In this paper we 
select the minimization of the energy spent by a swimmer to swim at fixed 
(given) average velocity c/T as the optimality criterion. 

2.1 Exterior Hydrodynamic Interactions 

Independently of the optimality criterion that one chooses, any swimmer 
needs to satisfy certain constraints which are required for the problem to be 
well posed and for it to actually represent a physically admissible swimmer. 

An organism is said to be swimming if it advances while its shape changes 
periodically in time by following the hydrodynamics of the surrounding fluid, 
in the absence of external forces or torques. 

More precisely, then, a low Reynolds number swimmer satisfies the fol- 
lowing key properties: 

• the velocity and force density at the surface of the swimmer are linked 
by the exterior Stokes equations satisfied by the surrounding fluid; 

• the shape of the swimmer changes periodically in time, driving the 
surrounding flow; 

• the force is generated internally by the swimmer, that is, the total 
external force and torque that act on the swimmer are zero. 

On the exterior part of the swimmer, this is expressed by the following 



4 



system of equations, satisfied for each t in [0, T]: 

-T]Au(x,t) + Vp(x,t) = -V-<r = in $l 3 \Q(t) (la) 

V-u(aj,0 = ° in^ 3 \fi(t) 

w(cc, £) = u(je, £) on 
lim |it(cc, £)| = 

| EC | — ¥OD 

lim |p(aj, £)| = 

\x\— >oo 



— = «(X,t) VXGT(t) (lb) 

T(T) = HT(0) + / [ vdTdt = RT(0) + ce 



an 

r(t) 



(lc) 



A ern = 

T(t) 

where u and p are the velocity and hydrodynamic pressure fields in the 
exterior domain Ji 3 \f2(£), r\ is the viscosity of the fluid, v is the velocity at 
the surface of the swimmer, X is a material point on the surface T(t) and e 
is the direction along which the swimmer is moving. 



The set of Equations (la) describes the conservation of momentum and 
mass in the Stokes fluid surrounding the swimmer, with zero boundary con- 
ditions at infinity. 



In Eq. (lb) we grouped the no-slip boundary condition (each material 
particle of the surface of the swimmer is following the evolution of the flow, 
i.e., no slip occurs between the fluid and the swimmers) and the periodicity 
of the swimmer shape, where c := f^v ■ edt is a shorthand notation for 
the distance swan by the swimmer. Notice how in the general case the 
configuration of the swimmer at time T is a rigid motion of the original 
configuration T(0), where R indicates a rotation and ce is the translational 
vector. 



Eq. (lc) states that forces and torques are generated internally by the 
swimmer. Here cr = r/(Vu + V T u) — pi is the stress in the fluid, with i" the 
identity. 
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In the axisymmetric case, the second of Equations (lc) is automatically 
satisfied, and we will no longer refer to it. Moreover, due to the symmetry 
of the system, no rotation is allowed and -R is the identity matrix, which we 
will omit from now on. Notice that this also implies that both velocities and 
generated force densities are axisymmetric. 

Knowledge of the velocity v(x,t) on the boundary is enough to ensure 
existence and uniqueness of a solution u(x,t) in the entire external domain, 
provided that T(t) is maintained regular enough throughout time: a bound- 
ary T(t) of class C 1 and boundary data v e Hz(r(t)) guarantee existence 
and uniqueness of the exterior solution u at time t, see, for example Ref. [?]. 

2.2 Global Hydrodynamic Interactions 



In Equations (la) we do not describe what happens inside the domain Q(t). 
While in reality micro swimmers have a rather complex internal structure, 
we will make two simplifying assumptions: 

• all movement capabilities of the organism are achieved through changes 
on the surface T(t) only; 

• the part of the swimmer enclosed by T(t) contains fluid identical to the 
one in which the swimmer is immersed, hence the same equations hold 
in the interior and in the exterior of Q(t). 

Under these assumptions, a swimmer is fully defined if we either assign 
a compatible movement of its boundary or by specifying the distribution of 
stress jumps that the swimmer is sustaining on T(t). 

The underlying principle we are exploiting here is that movement is 
achieved via hydrodynamic interaction between the surface of the swimmer 
and the surrounding fluid. Due to the complexity and diversity of the internal 
structure of biological swimmers, we don't attempt to model it explicitly. We 
do allow, however, for the presence of a fluid inside Q(t) because this is the 
typical case for bio-swimmers. We assume for simplicity that the viscosities 
of inner and outer fluids are the same because this allows one to solve the 
global Stokes problem using only the single layer potential. More precisely, 
the following relationship holds (see, for example, Ref. [8]): 

Snrjuix) = - J^S(x-y) i<r{y)n{y)\ dV(y), x e 9£ 3 (2) 
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where <S is the three-dimensional free-space Stokeslet 



5 ( r ) : = o + ^rnr ( 3 ) 



and the fluid flow is solved for in the entire space 3ft 3 . 

To simplify the notation a little, we introduce the definition of a Neu- 
mann to Dirichlet ADr and its inverse Dirichlet to Neumann map Z2Vr as the 
operators that, for each stress jump distribution / = [crn] on T(t), return 
the velocity distribution v generated on V that satisfy Eq. ([2]), and vice versa, 
i.e.: 

AD r / = v LN r v = f 

8n V v(x) = - j S{x- y)f(y) dT{y), xeT. (4) 

In the examples we consider, we always work with C 1 boundaries T. 
Therefore Eq. M) makes sense for / 6 H~^(T) and gives v G H^(T). 

We remark here that both the Dirichlet to Neumann map and the Neu- 
mann to Dirichlet map §4§ are translation invariant with respect to the con- 
figuration r(£), i.e., 

r 1 = r + d^AD ro / = AD ri /, LN To v = LN Tl v V«, V/, (5) 

for any constant translation vector d, i.e., the same stress jump applied to a 
translation of F produces identical velocity distributions on T itself. 

Given the above hypothesis, a biological swimmer can then be modeled 
entirely as a closed surface that is capable of transmitting forces to the ex- 
ternal (as well as to the internal) fluid in the form of stress jumps. The time 
evolution of the swimmer itself is then given by transporting the boundary 



with the flow field u(x, t) as in Eq. (lb), and the system of equations satisfied 
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by the unknown axisymmetric v and / becomes simply 
£N m v = f 



(6a) 



dX 

~dt 



v(X,t) 



r(T) = r(o) + I J v dr dt = r(o) + ce 



VX G T(t) (6b) 
(6c) 



/ = o, 



(6d) 



r(t) 



for t G [0, T], where Eq. (6a) effectively replaces the set of Equations ( la) and 
we already made the simplifications derived from the symmetry of the system: 
no external torque can be generated by axisymmetric force distributions, and 
the only possible motion is a translation along the axis of symmetry. 
More compactly, system ^ can be written symbolically as 



MS)-* 

r(r) = r(o) + ce . 



o 



Vt G [0,T] 



(7) 



Using cylindrical coordinates, the configuration of the swimmer, its force 
density distribution and its velocity field on T(t) can all be expressed as 
a rotation around the axis of symmetry (i.e., the x axis) of one-parameter 
(vector valued) functions (see Figure [l] for an example configuration). 

In particular we fix the domain T> of the parameter s to be the interval 
[0, 1]. The three-dimensional configuration V is the image of the function X, 
given by 

a(s)cos(9) \=R(6)X{s), 0g[O,2tt], sG[0,1], (8) 
<r(s)sin(8) J 

while the three-dimensional force density and the velocity fields are given by 



/(*,' 



/ Us) 

f*(s)cos(9) 
V A(s)sin(0) 



R(9)f(s), 0g[O,2tt], sg[0,1], (9) 
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Figure 1: An example of the instantaneous configuration of an axisymmet- 
ric swimmer on the (x, a)— plane. The three-dimensional configuration is 
obtained by a rotation around the a;- axis of the function X(s). 



and 

(u x (s) 
u a (s)cos(6) \=R(6)u(s), 6e[0,2ir], se[0,l], (10) 
u a (s) sin(0) 

where R(0) is defined as 

/I 

R(6) = cos(#) | . (11) 
\ sm(9) 

To avoid confusion, we will distinguish between three-dimensional func- 
tions u(x) and their axisymmetric counterparts u(s), identified with their 
section on the (x, a)— plane, by either using symbols with a tilde on top (e.g., 
S for the three-dimensional Stokeslet and S for the axisymmetric one) or by 
explicitly specifying the domain variable (x for three-dimensional functions 
or s for their (x, a) section). 
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Our approach is to assign T through a set of time dependent scalars 
£(t) G V C Sft^, so that for each time t, the curve X(£(£))(s) generates an 
admissible T. We focus on cases where the boundary T depends smoothly 
(e.g., analytically) on the parameters £, and remains of class C 1 for all time, 
X(s, t) is a non self-intersecting, C 1 curve on the (x, a)— upper half plane (a > 
0), and the surface obtained by its rotation around the x axis is the boundary 
r(t) of a domain Q(t) such that its complement 3ft 3 \ Q(t) is connected. 

The admissible shapes we treat in this paper are collections of simple 
closed curves in the (x, er)— upper half plane, or open curves in the same half 
plane which start and end vertically on the x axis. 



2.3 Optimal Swimming as a Constrained Minimization 
Problem 

While system ^ describes the relationship between the forces / that the 
swimmer applies to the surrounding flow and its consequent evolution v, it 
gives no information about the optimality of a given swimming stroke. 

A natural optimality criterion consists on minimizing the energy dissi- 
pated while swimming at a given average velocity c/T. In this sense the 
optimal stroke is the one that minimizes 

£{v):= f [ v-LNrvdTdt, (12) 
Jo Jv{t) 

where v satisfies the set of Equations ([6]), subject to the constraint 

F e (v):= [ [ v-edTdt-c = 0. (13) 
Jo Jv(t) 



In order to embed in the minimization problem also Eq. (6d), it is con- 
venient to express the entire problem in terms of shape changes rather than 
absolute velocities. 

We assume that the changes in the configuration T(t) of the swimmer 
happen only through a set of iV + 1 scalar functions of time which we identify 
as the shape variables i — 1 . . . N, and the location (p (t) on the axis of 
symmetry of a distinguished point. 

Moreover, the configuration T (£,<£>) is known for any admissible shape £ 
and displacement ip G 3ft and has the form 

r(£, ip) := r (£) + e<p, (£, (p) E 3ft w x 3ft, (14) 
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as exemplified in Figure [2] 

In particular, we assume that the configuration F is the image of a func- 
tion X of the parameter s6D, which identifies the location of each material 
point of T, i.e., 

X(£,<p,s) :=X(^s) + e<p, r(£, <p) = {X(£, <p, s), s G V} . (15) 




Figure 2: Example of configuration T(£,(p). 



Using Equations (6b) and (15), the velocity of the swimmer can then be 
expressed in terms of rate of shape changes £ and in terms of translational 
velocity (p as 

v{X(£(t),(p(t),s)) = — &(t) + (p(t) 



dX(£,s) 
<9& 



dip 



£i(t) + e<p(t) 



(16) 



--Ui(s)Ut) + e(p(t), Vs in V, Vt G [0,T], 



where summation is implied on repeated indices. Here e is a unit vector di- 
rected along the axis of symmetry, while Ui is the change in the configuration 
r(t) at the position X(£, ip, s) due to a variation of the shape parameter £j. 
Using this formalism, the periodicity condition (6c) becomes simply 

£(0)=f(T). 
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Eq. (6d) implies that 

LN r v-edT = (17) 



r 



which can be rewritten using Eq. (16) as 



J LNriuii + ecp)-e = • £ + M(0 <f = 0, (18) 



where 



Ni(0 := J I^vUi- edr 
M(f) := / HVre-ecir. 



(19) 



r 



It is now possible to write ip as a function of £ 



which allows us to express the energy dissipated by a given stroke £ as 

£(£):= / T G(0£-£dt, (21) 



where 



(0 :=y LNrWi-Wjdr, w t := Ui + V^e. (22) 

A natural space where we can look for solutions £ of our problem is then 
given by the Sobolev space 

V := {£ G H\0,T) N , s.t. f(0) = f(T)} (23) 

of periodic i/ 1 functions of N components. Optimal swimming is then given 
by the solution £ G V of 

mm [ T G(^-Ut =:£(£) (24a) 
? ev Jo 



T 







subject to / V(f) • £ dt - c =: J" c (£) = 0. (24b) 
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In principle, additional constraints could be added to Eq. (24b) to fix, 
for example, that the swimmer initial and final shape are given, or that 
the shape never reaches out of given physical bounds (due, for example, to 
construction constraints in artificial swimmers or physiological constraints in 
biological swimmers): 



£(0) = £(T) = £ 

a < m < & 



Wt e [o,T]. 



(25a) 
(25b) 



In Ref. [5] and Ref. [6], problem (24) with the addition of constraint ( |25a ) 
has been considered for some special model swimmers. In this case, the space 
of admissible shapes is 

V := {£ G H\0, T) N , s.t. £(0) = £(T) = £ } (26) 

and the following hypothesis guarantees the existence of a solution to prob- 
lem pi}: 



There exists £ G Vo which satisfies (246) and such that £ (£) is finite, (27) 
as shown below. 

Proposition 1 Assume G and V are continuous functions and that there 
exist positive constants a, such that 



3a > 0, (3 > 0, such that , P\r)\ 2 > r] ■ G(£)t? > a\r]\ 2 Vr] 6 



(28) 



uniformly for admissible shapes £ G Vo- Then there exists a solution to 
problem (24)- 

Proof 1 The proof is a straightforward application of the Direct Method of 
the Calculus of Variations. Indeed, under the given hypotheses, the func- 
tional to minimize is coercive and lower semicontinuous with respect to weak 
convergence in (H l (0,T)) N . Moreover the constraint is continuous with re- 
spect to this topology. Condition (21) ensures that the set over which the 
minimization is done is not empty. 



Remark 1 Condition (21) is nontrivial. It requires that the swimmer can 
indeed cover the prescribed distance, using finite energy. As explained in 
Ref. [5] and Ref. f^, this is a controllability condition equivalent to asking 
that the constraint (24b) is non-holonomic. 
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Remark 2 The uniform coercivity condition (28) is a very reasonable con- 
dition although very difficult to prove in a framework as general as the one 
presented here. What is however clear in view of its definition, is that for all 
£ G Vo the matrix <&(£) is positive definite. 



We will not pursue further the analysis of problem (24), which requires 



analytical verification of conditions such as (27) and (28) on a case- by-case 
basis. In the rest of our work, we will focus instead on strategies to find 



solutions to the optimization problem (24) numerically. 



3 Numerical Approximation 

The two main ingredients of the computer simulation of the swimmer prob- 
lem are an efficient numerical implementation of the Dirichlet to Neumann 
map Q and a strategy for the search of minima of £(£) in Eq. (24) satisfying 
the constraints. 

We solve the first problem using the boundary element method (with a 
custom code written in C++, using the deal . II library [HI 10] and based on 
the formulation presented in Ref. [8]). We address the solution to the con- 
strained minimization problem using the reduced space successive quadratic 
programming strategy (rSQP)[TT] provided with the Moocho package of the 
Trilinos C++ library. [ 



3.1 Axisymmetric Boundary Integral Equations for Stokes 
Flow 

Our boundary integral formulation for axisymmetric Stokes flow follows closely Ref. [8] : 
if we consider Stokes equation in free-space associated with a forcing term 
due to a Dirac mass centered in y and weighted with the force vector b, i.e., 

- r)Au(x) + Vp{x) = -V • cr = bS(x - y) in -ft 3 

V ■ u(x) =0 in 
lim \u(x)\ = (29) 

\x\— >oo 

lim \p(x) \ = 0, 

\x\— >oo 
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we can express the solutions u, p and cr using the free-space Green's functions 
<S, P and T: 



u(x) = Six — y)b 

V ' 87TT] V y; 

= irP( x -y) b 

crfa?) = — fix - y)b, 



(30) 



where 



5 A 
\r 



]l_ _|_ r » r i 



T, 



ijk \ 



I' 3 

riTjTk 



(31) 



The axisymmetric approximation allows us to reduce the swimmer prob- 
lem to a boundary integral equation on a one-dimensional curve, i.e., the 
swimmer is allowed to generate only axial-symmetric surface forces, which 
are then transmitted to the fluid. 

We introduce the following notation for integrals on the surface T, iden- 
tified by the configuration X(s) = (X(s), cr(s)), s G [0, 1]: 



g dr = 2tt 



g(s)J x (s)a(s) ds 
g(s)ds, 



(32) 



X(s) 



where we indicated with Jx(s) 



— 
ds 



and the last line serves only to 



ease the notation. 

To solve Stokes equations in the axisymmetric case, we integrate analyt- 
ically the singular integral Eq. ^ around the axis of rotation, while keeping 
the point x on the half-plane 6 = 0. 

The azimuthal integration gives rise to a single layer potential kernel, 
which can be expressed using complete elliptic integrals of the first and sec- 
ond kind and produces a linearly coupled set of one-dimensional integral 
equations. 
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We can write the resulting equation as 



8irr]V a (q) 



S a/3 (X(q),X{s))fJs)ds, 



(33) 



X(s) 



where the indices a and are either x or a, to indicate the distance along 
the axis of symmetry, or the distance from the axis itself. 

The axisymmetric kernel «S represents the single layer integrals of a ring 
of singularities passing through the point X(s, 0) = R(0)X(s): 



S(X,Y) 



2 " r S xx (S yx cos 9 + S zx sin 9) 
<Sx V (S yy cos 9 + S zy sin 6) 



(34) 



where S = S(R(0)X, R(9)Y) is the three-dimensional Stokeslet given in (31 ), 
evaluated at the points R(0)X and R{6)Y . 

The above integration can be expressed in terms of complete elliptic in- 
tegrals of the first and second kind, which are defined, respectively, as 



K(x) :-- 



tt/2 



d# 



'l — x 2 sin 2 6 
which is a logarithmically singular integral, and 

f7r/2 



da, x G [0, 1), 



(35) 



r 1 i 

E{x):= VI -x 2 sin 2 fldfl, x G [0, 1], 

Jo 



(36) 



which is a bounded integral. 
Indeed, if we define 



00 



r := a/(x - x ) 2 + (a - a f 



Ax = x — xq, 



and 



4(T(Tn 



fx - x ) 2 + {a + cr ) 
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then the axisymmetric Stokeslet is given by (see Ref. [8]) 









= 2k J — 




V °o 




kAx 








a \ 




kAx 






O"0 \ 



K(k) + Ax^ 




K(k) - (o% - a 2 + Ax 2 




K(k) + (ol -a'- Ax 2 )^) 



a a y cr Q 




(a 2 + a 2 + 2Ax 2 )if(fc)- 



(2Ax 4 + 3A;r 2 (cx 2 + a 2 ) + (a 2 - a 2 ) 2 



E{k) 



(37) 



where the diagonal entries are singular and behave like S xx ~ <S CT<7 ~ — 2 ln(r) 



3.2 Spatial Discretization 

We employ an iso-parametric finite dimensional representation for the spatial 
discretization of the configuration X, and of the velocity and force density u 
and /. We call Vh = span{r l }i<j<M the finite dimensional space that we use 
for the discretization. To ensure regularity and smoothness of the boundary 
T for any t, we choose a discretization based on cubic bell splines. 

In particular at any given time we identify elements of Vh with the coef- 
ficients of M dimensional vectors, as in 

X h (s) = XW\s) 

V h {s)=V i T i (s), 

where we used the summation convention, and i goes from 1 to M. The su- 
perscripts here indicate the indices of the basis functions that span the space 
Vh- When confusion cannot arise, we will use X h without "(s)" to indicate 
the 9ft M vector of components X\ i = 1 . . . M which uniquely identify the 



vector valued function Xh(s) in Vh as in (38). 

To take into account collections of shapes, it is possible to introduce con- 
trolled discontinuities in the basis functions r*(s), by overlapping a sufficient 
number of nodes in the knot-span that generates the B-Spline basis, as ex- 
plained, for example, in Ref. [T3] . 
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The Galerkin approach to the boundary element method can be under- 
stood by multiplying Eq. (33) on both sides with test functions (or virtual 
forces) r and integrating on T 




v(s) ■ t(s) ds = 



~[ r{q)~f S(X(q),X(s))f(s)dsdq 
Jx(g) 8nr] J x{s) 

Vr e C°°([0,1]) 2 , (39) 

which, on a discrete level, yields the discrete version of the Dirichlet to Neu- 
mann map: 

f h = A^Muh, (40) 
where both matrices A and Ai depend on Xh(s) and are defined as 

A ij = - [ r\q) ■ / S(X h (q),X h (s))T>(s)dsdq, (41) 
Jx h (q) onr] J Xh{s) 

and 

M ij = [ T j (s) ■T i (s)ds. (42) 

Numerical computation of the matrix entries of both A and Al is per- 
formed using high order Gauss quadrature formulae, as well as Gauss-Log 
quadrature formulae for the diagonal entries of A, in order to properly take 
into account the logarithmic singularity of the axisymmetric Stokeslet (see, 
for example, Ref. El H5]). 

Notice that, once we have two vectors f h and Vh in !ft M that represent 
objects of Vh, we can compute their integrals on Xh(s) by scalar product in 
sjjm through the matrix At: 

/ f h (s) ■ v h (s) ds = f h ■ Mv h = v h ■ Mf h = v { M ij f, (43) 

Jx(s) 

where in the left hand side (■) has the meaning of scalar product between 
vector functions of 3ft 2 , while on the right hand sides (•) should be read as 
the scalar product in 3? M . 
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Once we know how to discretize a domain T and write a discrete version 
of the Dirichlet to Neumann Map for any configuration, we simply define the 
discrete shape basis functions as the counter part of dX(£,(p)/d£, i.e., 



% ' (44) 



u\(s) = ' 1 = 1. ..n 



e h (s) = e{s), 

where X(£, s) indicates a point on the (x, a)— half-plane of the configuration 

r (0- 

The discrete shape basis functions u\(s) are simply constructed as the 
b-spline interpolation of their continuous counter parts, which are a datum 
of the problem, while the discrete basis function e^(s) for the displacement 
vector space is simply equal to the discrete representation of the constant 
unit vector directed along the axis of symmetry. 

Given N + 1 arbitrary time functions (£(£), f{t)), it is then possible to 
build the discrete configuration Xh(£(t),(p(t)), the discrete basis functions 
for the rate of shape change u\ and to compute the Dirichlet to Neumann 
map for arbitrary data defined on the discrete curve. 

In particular we can compute the discrete versions of V and G as 

■= " (A->Me h )-Me h ' (45) 

and 

G ij>h {0 := (A-'MiuUO +V ilA (£)e A )) • M(u{(£) +V j j i (Z)e h ). (46) 
The semi-discrete formulation of the optimal swimming problem is simply 



obtained by substituting in Eq. (24) V and G with their spatial discrete 
counterparts: 

mm [ T G h {$i-idt=:£ h {0 (47a) 
« ev Jo 

subject to / V h (£)-idt-c=:T c h (O = 0. (47b) 
Jo 
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3.3 Time Discretization 

The dependency of the domain T (or, equivalently, of its axisymmetric rep- 
resentation X) on the time variable t is only through the shape and position 
functions ^(t) and (pit). In particular we rewrote the problem in a way that 
makes it only dependent on the shape parameters £ G V. 

We approximate also these functions with cubic b-spline functions, with 
a possibly different dimension Q. In particular we have 

V h = G (span^t)}^)^ , s.t. £,(0) = &(T)} , (48) 

where we use Greek letters to distinguish basis functions referring to the time 
variable t from those referring to the space variable s and again we can write 



Zh(t) i = &'Y a (t) i = l...N, (49) 

where summation is implied on the repeated indices a which go from 1 to Q. 

Notice that a numerical approximation ^ is now a vector of dt Nx ®, which 
we identify either with two indices £f , or with a single capital letter index 
I — 1 ... N x Q, where 

I:=i + (a- m & ■= {{£f }f=i}L = fe}f=i Q , (50) 

that is 



£l> ^2> • • • ) £iV; £l ) ^2) • • • ) £iV! •••) £l > £2 ; • • • ; £iV • (51) 



3.4 Discrete Optimal Swimming 

With the discretization Vh of the infinite dimensional space V, we can reduce 
the minimization problem to a finite dimensional one defined on §t® xN , where 
Q is the total number of basis splines selected for the discretization of the N 
shape parameters in time, and the fully discrete problem becomes 

min £ h (S h ):=<S(Sh)Sh-Sh (52a) 
subject to J- h c (a) := • Ch ~ c = 0, (52b) 
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where the matrix <£(&) in sr(XxQ)x(NxQ) and the vector 3^) in sr(NxQ) are 
defined as 

T G ij (utm a (t)-i f) (t)dt 

(53) 



and the relationship between the indices I, J and i,j,a,(3 is given by (50) 
and J = j + (/3- 1)N. 

The problem is now written in the framework of classical finite dimen- 
sional optimization. Indeed, if we define the nabla operator V applied to a 
function of the vector variable ^ as 



(VFfo))/ :- 



(54) 



we can then write the Lagrangian C(C,h, A) and its derivatives associated with 
problem (52) as 



£(&,A) 
V£(&,A) 
V 2 £(6>, A) 



-S h (Ch) + A^ c (a) 

=v 2 ^(e,) + Av 2 ^(a). 



(55) 



The first and second-order necessary Karush-Kuhn- Tucker (KKT) opti- 
mality conditions for a solution (£/,, A) to (52) are: 



V£(a, A) = V5 h (&) + AV7£(&) 

m h ) = 

V 2 £(^ h ,X) v - v >0 V»ieV f 







(56) 



also known as linear dependence of gradients, feasibility and curvature con- 
ditions. 

A popular class of methods for solving non linear constrained minimiza- 
tion problems is successive quadratic programming (SQP) (see, e.g., Ref. [IT]), 
which is equivalent to applying Newton's method to solve the the minimality 
conditions (56). 



At each Newton iteration k for (56), the linear subproblem (also known 
as the KKT system) takes the form 

V£( fc ) 



W A 
A T 



d 



(k) 



(57) 
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where 



d 


._ Ak+l) Ak) 


d\ 


■- _ \(k) 


W 


:=V 2 £(ef,A«) 


A 


:= V^ c (£f ) 




A«) 


7T(fc) 
J c 





(58) 



and we use superscripts between parenthesis (k) to refer to the k-th Newton 
iteration, not to be confused with the index of the b-spline basis functions. 
Once we obtained a new estimate of the solution 1 , A < - /c+1 - ) ), the error 



in the optimality conditions (56) is checked. If these KKT errors are within 
some specified tolerance, the algorithm is terminated with the optimal solu- 
tion. 

If the KKT error is too large, the functions and gradients are then com- 
puted at the new point and another KKT subproblem (57) is solved 



which generates another increment d, until convergence or failure. 

A successful application of the SQP method requires one to provide exact 
information about the hessian of the Lagrangian W. While this ensures 
second order convergence of the SQP method, it is not feasible in many 
practical cases, including ours, where the direct computation of W is overly 
expensive. 

To address this issue, we use a reduced-space SQP (rSQP) method. In 



rSQP methods, the full-space QP subproblem (57) is decomposed into two 
smaller subproblems, where the optimization vector variable ^ is split into 
dependent and independent variables, by using the linearized equality con- 
straint as a mean to express the dependent variables (often referred to as the 
control variables) in terms of the independent ones (also referred to as the 
state variables). 

Using this decomposition, the two subproblems can be solved effectively 
by using an approximation of W. We use the implementation of the rSQP 
method included in the Moocho package of the Trilinos library. 
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4 Numerical Results 



The goal of this section is to present some results obtained with our code. In 
section 4.1 we present our results for the optimal stroke of the Golestanian 
swimmer. These should be confronted with the results already presented 
in Ref. [5] and Ref. [6] with a completely different approach both for the 
resolution of Stokes equations and for the minimization of the expended 
energy, which we take as our reference results. 

While from the quantitative point of view we observe some marginal 
differences, qualitatively the method presented in this paper is capable of 
reproducing the optimal strokes obtained in the previously mentioned works. 

We should mention here that the main differences lay in the fact that in 
those works the approximation of the V(£) and G(£) fields was performed 
using a Finite Element Method for the simulation of axisymmetric Stokes 
Flow in a box that was taken large compared to the dimension of the swim- 
mer, but of finite dimension. In comparison, the method we present here is 
based on the use of Boundary Integral Approximation, which should yield 
more accurate results on infinite domains such as the one we are interested 
on. 

As a second remark, we observe that the optimal stroke in Ref. [3] and 
Ref. [6] is obtained by writing explicitly the Euler-Lagrange equation for 
the constrained optimization problem, and by solving the resulting system 
of ODEs with an explicit Runge Kutta method, coupled with a shooting 
approach in order to enforce the periodicity of the strokes. 

In this paper, on the other hand, the solution of the minimization problem 
itself is left to an external library which explores the space of shapes in 
a very efficient and flexible way, allowing us to solve problems that could 
not be addressed before. In particular we can solve problems in which we 
include inequality constraints on the shape variables £, and we can release 
the constraints on the initial and final shape, letting the minimizer find the 
optimal starting shape for our swimmers. 

In the experiments that follow, we use water at room temperature (20° 
C) as the surrounding fluid, and we express lengths in millimeters (mm), 
time in seconds (s) and weight in milligrams (mg). Using these units, the 
viscosity of water is approximately one (ImPas = lmg mm' 1 s' 1 ), and the 
energy is expressed in pico- Joules (lpJ = lmg mm 2 s~ 2 = 10~ 12 J). 

The tests were performed on a Macbook Pro, with 2.16 GHz Intel Core 
2 Duo processor and 2 GB 667 MHz DDR2 SDRAM. The average running 
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Figure 3: Definition of the Three Sphere Swimmer. 



time for the experiments presented in Sec. 4.1 is about 30 seconds, with a 



maximum memory consumption of 64MB, while the experiments of Sec. 4.2 



required on average 10 seconds, with roughly the same memory consumption. 

The low running time is possible thanks to the use of cubic B- Spline, 
which allows one to obtain accurate solutions using very few degrees of free- 
dom. In all the experiments we present, we used 15 control points on each 
portion of the swimmer for the spatial discretization (i.e., 45 total control 
points for the three sphere swimmer and 30 for the stick and donut swimmer) 
and 15 control points for the time discretization of each shape and position 
variable. 



4.1 Three Sphere Swimmer 

The three sphere swimmer is among the simplest axisymmetric swimmers 
one can think of. It consists of three linked spheres that can only vary their 
reciprocal distance. The swimmer is completely defined once we have the 
radius of the three spheres (which we assume to be the same for simplic- 
ity), and the location of the three centers of the spheres on the axis of the 
movement (one size parameter and three positional variables). 

Equivalently, we can use the more convenient representation given by the 
location of the center of the central sphere (ip) and the distances between 
the two lateral spheres and the central one (£i and £ 2 )> as shown in Figure [3] 
These two different representations are clearly equivalent (see, e.g., Ref. [5]). 
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Using this convention X(£, <p, s) can be parametrized as 

*«^K"^-(3t:r s)> ) ifse|o - i/3) 

which implies that the basis functions for shape changes Ui := dX/d^i are 
defined as 

ui(s) = -X[o,i/ 3 )(s)e, u 2 (s) = X(2/3,i](s)e, (60) 

where Xa{s) is the function which is one if s G A, and zero otherwise. We 
constrain the shape variables & to be in the interval [0,6a]. Notice that 
in this paper £ is the touching distance between the spheres, and not the 
distance between the centers of the spheres, as in Ref. [S]. 

In Figure [4] we show the corners of the box [0, 6a] 2 in which the three 
sphere swimmer shape is constrained to stay, for a swimmer of radius a = 
.05mm. 

m — e e — e — • 

00 

Figure 4: Extremal shape configurations for a Three-Sphere swimmer of 
radius a = .05mm, from left to right and from bottom to top: £ = (0,0), 
£ = (6a, 0), £ = (0, 6a) and £ = (6a, 6a). 

The basis function for the change in position, dX/dip is simply equal 
to the unit vector e. In Figure [5j we show the three basis functions that 
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allow one to fully describe the system, together with their Dirichlet to Neu- 
mann maps, evaluated at the configuration £ = (.05mm, .05mm), which is 
equivalent to asking that the distance between the centers of the spheres is 
3a. 




Figure 5: Basis function Uj(s) and e (left) and their Dirichlet to Neumann 
maps (right) for £ = (a, a). 




Figure 6: The force-free basis function Wi(s) (left) and their Dirichlet to 
Neumann maps (right) for £ = (a, a). 

Once we have the various basis functions, it is easy to combine them 
linearly and obtain a basis for force-free movement. This is what Figure [6] 
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shows, where the basis functions Wi = Ui + V,e are plotted, together with 
their Dirichlet to Neumann map. Notice that, by construction, the integral 
of the forces on the configuration V yields zero to machine precision. 

Minimizing the expended energy and forcing the displacement to be a 
given datum, we obtain the same qualitative results as in Ref . [5] and Ref . [B] . 
The approach we propose in this paper, however, allows us to study also the 
problem of finding the optimal stroke without assigning an initial (and final) 
shape. 

The closed paths in shape space for target displacements of .01mm and 
.001mm can be seen in Figure [7] for both the case where the initial shape is 
fixed to be (.2mm, .2mm) and for the case where no constraints are imposed. 

— i 1 1 1 1 

c = 0.01mm 

c = 0.001mm 

O ; 

i i i i i 

0.1 0.2 0.3 0.1 0.2 0.3 



Figure 7: Path in shape space describing optimal strokes for three sphere 
swimmer of radius .05mm, swimming .01mm and .001mm in 1 second, im- 
posing the initial shape £ = (.3mm, .3mm) on the left, and without imposing 
an initial shape on the right. 

In Figure [8] we can observe how the optimal stroke can be separated into 
a propulsive part and a recovery part, where the propulsive part is the one 
where the velocity is bigger and positive, while the recovery part is the one 
where the velocity is smaller and negative. 

We would like to emphasize here that, in general, if we fix the initial and 
final shapes, the velocity is free to be discontinuous at t — kT with k integer, 
however the truly optimal stroke (the one selected without impositions on 
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c = 0.01mm 



c = 0.01mm 




Figure 8: Diagram of translational velocities for three sphere swimmer of 
radius .05mm swimming .01mm in 1 second imposing the initial shape £ = 
(.3mm, .3mm) on the left, and without imposing an initial shape on the right. 

the starting shape) is one where also the velocity is continuous, as shown in 
Figure [7] on the right. 

4.2 Stick and Donut swimmer 

We now present a new model swimmer, which simulates the swimming mech- 
anism of many biological organisms made of a body and a propulsive appara- 
tus consisting of appendages that change shape to induce propulsion. A nice 
example is Chlamydomonas Reinhardtii, a unicellular eukaryotic alga with a 
body size of roughly 10 /im, that swims by executing a movement with its 
two flexible flagella which is closely reminiscent of the breast stroke of an 
Olympic swimmer. An axisymmetric version of this swimming style, though 
at larger scales, is that of the jelly fish. In the model swimmer we propose, 
the body is schematized with a cylinder capped with two half-spheres at the 
ends (the Stick) and the propulsive apparatus is schematized with a torus 
with elliptical cross-section of variable major and minor radii (the Donut). 

Figure [9] shows a section of the swimmer on the (x, <j)-half plane, for 
6 = 0. The radius of the stick is set to be Ro, and all other dimensions are 
scaled with respect to this one, so that the touching distance between the 
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stick and the donut is fixed to 3R /2 and the volumes of the stick and the 
donut are fixed and equal. 




Figure 9: Definition of the "Stick and Donut" swimmer. 

Swimming is achieved by translating the center of the donut along the 
direction of the stick (non-dimensional shape variable £ 1; constrained in the 
interval [0, 1]), so that the center of the donut is always inside the "body" of 
the swimmer), and by varying the radii of the donut section. 

For the non-dimensional shape parameter £2 = 0, the horizontal radius 
Ri of the donut section is equal to 5R /2, i.e., it has the same length of the 
stick without the half spheres, while for £ 2 = 1, the horizontal radius Ri is 
equal to Rq/5. 

The vertical radius of the swimmer and the vertical position of the center 
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of the ellipse section are adjusted automatically as a function of £ 2 in order 
to maintain the volume of the donut constantly equal to the volume of the 
stick and to maintain the distance between the donut and the stick constant. 

The absolute position of the swimmer on the axis of symmetry x, is given 
by the non-dimensional variable (p, so that when (p is equal to one, then the 
swimmer has moved of one entire body length (7Rq). 







r 



.„//lt'ti\u... 



Figure 10: Basis functions ti;(s) and e(s) (left) and their Dirichlet to Neu- 
mann map fi(s) and f e (rigth) for £ = (.5, .5). 



Figures 10 and 11 show the basis functions it, and IR^e and the corre- 



sponding force-free basis functions Wi with their Dirichlet to Neumann maps. 



From Figure [TT] it is evident how changes in the shape parameter £2 do not 
induce much displacement along the axis of symmetry. 

We constrain the swimmer shape variables to be included in the non- 



dimensional box [0, 1] , whose corners are shown in Figure 12 A square 



stroke that explores in a clockwise manner this shape space is presented in 



30 




Figure 11: Force free basis functions Wi(s) (left) and their Dirichlet to Neu- 
mann map (rigth) for £ = (.5, .5). 



Figure 13 



The fact that changes in the shape parameter £2 have little effect on the 
overall displacement of the swimmer can also be inferred from the left and 



right sides of the square stroke, in the right part of Figure 13, which shows 
how the velocity ip of the swimmer due to changes in £ 2 are negligible when 
compared with the velocity due to changes in £1 (top and bottom sides of 
the square stroke). 

The separation of stroke cycles into a power or propulsive phase, in which 
appendages have maximal hydrodynamical resistance, and a recovery phase, 
in which the swimmer tries to minimize viscous drag forces on its propulsive 
appendages, are very common in nature. A typical example, schematically 



depicted in Figure 14, is the ciliary stroke cycle. 



Our model swimmer mimics this behavior by alternating propulsive phases 



(top side of the square stroke in Figure 13, i.e., £1 varying from zero to one 



while £2 is close to one) with a recovery phases (bottom side of the square 



stroke in Figure 13, i.e., £1 varying from one to zero while £2 is close to zero). 

In order to compare the stick and donut swimmer with the three sphere 
swimmer, we set the radius Rq of the swimmer such that the total volume 
of the stick and donut swimmer is equal to the volume of the three sphere 
swimmer with radius .05. This gives a radius Rq of approximately .034mm. 
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Figure 12: Extremal shape configurations for the stick and donut swimmer: 
from left to right and from bottom to top: £ = (0,0), £ = (1,0), £ = (0, 1) 
and f = (1,1). 



Figure 15 shows two optimal strokes found with our method which lead to 
displacements in Is of .01mm and .001mm. In the first case, the energy con- 
sumption is around .12QpJ, while in the second case it is .OlOpJ. This should 
be compared with the energy consumed by an "equivalent" three sphere 
swimmer, namely, a swimmer with the same volume as this one, swimming 
at the same average velocity. The three sphere swimmer energy consump- 
tion to swim .01mm in one second is .183pJ, while its energy consumption 
to swim .001mm in one second is .018pJ. 

In the long distance, the stick and donut swimmer is about 45% more effi- 
cient than the optimal three sphere swimmer, while for the shorter distance, 
the increase in efficiency is about 75%. 

A collection of animations referring to the optimal strokes of the swim- 
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Figure 13: Path in shape space describing a square strokes for stick and 
donut swimmer with base radius of .034mm, swimming .034mm in 1 second 
(left), with its propulsion diagram (right). 



mers presented in this work (Figures [F] and [IB} can be viewed on-line on the 
home page of the corresponding author. jTTj 
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